
#若在此单独运行改程序需设定框内参数
##################################################################################
# city_code<-'qd'   #目标城市
# pty<-11           #产品类型
# minmum_city<-30   #城市最小样本量设定
# minmum_dist<-30   #行政区最小样本量设定
# date_base<-24133  #基期设定为2011年1月
##################################################################################

city_list<-read.csv('OLS/city_list.csv')
city_name<-city_list%>%filter(city_cd==city_code)%>%select(one_of('city_name'))%>%unique
city_name<-city_name[[1]]%>%as.character


#抽取数据
source('R/load_spatial_data.R')

#ha_price与ha_phase数据处理
library(dplyr)
price_rent<-subset(tabln.vec$ha_price,proptype==pty&ymid>=date_begin&ymid<=date_end)[,c('ymid','ha_code','rentprice','rentbldgarea')]%>%na.omit
# phase<-tabln.vec$ha_phase%>%
#     filter(!is.na(buildyear))%>%
#     group_by(ha_code)%>%
#     summarise(buildyear=round(mean(buildyear)))%>%as.data.frame

#关联POI与info数据
ppi_rent<-merge(price_rent,ha_info.sp,by='ha_code')

#城市级别模型建立
cat(city_code,'model building...','\n')
myvar_city<-names(ppi_rent) %in% c('city_code','ha_code','proptype','rentcount','name','ha_cl_code','ha_cl_name','dist_code','dist_name','x','y','volume_rate','greening_rate')
rent<-ppi_rent[!myvar_city]

#城市最小样本量设定
tab.date<-table(ppi_rent$ymid)%>%as.data.frame%>%subset(Freq>=minmum_city)

date<-sort(unique(tab.date$Var1))%>%as.character%>%as.numeric
len<-length(date)

if(len==0){
    len_blank<-len_blank+1
    city_blank[len_blank,]<-city_name
}
if(len==0) next  #如果每月样本量小于30，跳出本次循环

fit.vec_city<-vector(mode='list',len)
names(fit.vec_city)<-date
for (i in 1:len){
    rent.new<-subset(rent,ymid==date[i])
    fit.vec_city[[i]]<-lm(log(rentprice)~.,data = rent.new[,names(rent.new)!='ymid'])
}

#计算城市HPI（拉氏指数）
cat('compute HPI of',city_code,'\n')
data_base<-rent[1,-1]

#若为新城市则用最早6个月样本量平均特征作为标准建筑
if(min(rent$ymid)<=date_base_end){
    data_base1<-subset(rent,ymid>=date_base_begin&ymid<=date_base_end) #基期设定
}else{data_base1<-subset(rent,ymid<=min(rent$ymid)+5)}

len_base<-length(names(data_base))
for(ii in 1:len_base){
    name_base<-names(data_base)[ii]
    data_base[,name_base]<-mean(data_base1[,name_base],na.rm=T)
}

result<-as.data.frame(cbind(dt<-rep(0,len),price.pre<-rep(0,len),r2<-rep(0,len),FBPI<-rep(0,len)))
names(result)<-c('DT','PRE','AdjR2','FBPI')
for(j in 1:len){
    price.pre<-predict(fit.vec_city[[j]],newdata=data_base)
    fit.sum<-summary(fit.vec_city[[j]])
    fit.squ<-fit.sum$adj.r.squared
    result[j,1]<-date[j]
    result[j,2]<-exp(price.pre[[1]])
    result[j,3]<-fit.squ
    result[j,4]<-result[j,2]/price.pre_base*pts_base
}

#行政区级别模型建立
myvar_dist<-names(ppi_rent) %in% c('city_code','ha_code','proptype','rentcount','name','ha_cl_code','ha_cl_name','dist_code','x','y','volume_rate','greening_rate')
ppi_rent1<-ppi_rent[!myvar_dist]
dist.name<-as.vector(unique(ppi_rent1$dist_name))
len_dist<-length(dist.name)
fit.vec_dist<-vector(mode='list',len_dist)
names(fit.vec_dist)<-dist.name
for(m in 1:len_dist){
    dist_cd<-dist.name[m]
    cat(dist_cd,'model building...','\n')
    rent_dist<-subset(ppi_rent1,dist_name==dist_cd)
    
    #行政区最小样本量设定
    tab.date_dist<-table(rent_dist$ymid)%>%as.data.frame%>%subset(Freq>=minmum_dist)
    tab.date_dist1<-tab.date_dist$Var1%>%as.character%>%as.numeric
    date_dist<-sort(unique(tab.date_dist1))
    len_date_dist<-length(date_dist)
    
    if(len_date_dist==0) next  #如果每月样本量小于30，跳出本次循环
    
    fit.vec_dist1<-vector(mode='list',len_date_dist)
    names(fit.vec_dist1)<-date_dist
    for (n in 1:len_date_dist){
        rent_dist.new<-subset(rent_dist,ymid==date_dist[n])
        fit.vec_dist1[[n]]<-lm(log(rentprice)~.,data = rent_dist.new[,names(rent_dist.new)!='ymid' & names(rent_dist.new)!='dist_name'])
    }
    fit.vec_dist[[m]]<-fit.vec_dist1
}

#计算行政区HPI（拉氏指数）
result.vec_dist<-vector(mode='list',len_dist)
names(result.vec_dist)<-dist.name
for(pp in 1:len_dist){
    dist_cd<-dist.name[pp]
    cat('compute HPI of',dist.name[pp],'\n')
    rent_dist<-subset(ppi_rent1,dist_name==dist_cd)
    
    #行政区最小样本量设定
    tab.date_dist<-table(rent_dist$ymid)%>%as.data.frame%>%subset(Freq>=minmum_dist)
    tab.date_dist1<-tab.date_dist$Var1%>%as.character%>%as.numeric
    date_dist<-sort(unique(tab.date_dist1))
    len_date_dist<-length(date_dist)
    
    if(len_date_dist==0) next  #如果每月样本量小于30，跳出本次循环
    
    data_base_dist<-rent_dist[1,names(rent_dist)!='dist_name'&names(rent_dist)!='ymid']
    
    #若为新行政区则用最早6个月样本量平均特征作为标准建筑
    if(min(rent_dist$ymid)<=date_base_end){
        data_base_dist1<-subset(rent_dist,ymid>=date_base_begin&ymid<=date_base_end) #基期设定
    }else{data_base_dist1<-subset(rent_dist,ymid<=min(rent_dist$ymid)+5)}
    
    len_base_dist<-length(names(data_base_dist))
    for(qq in 1:len_base_dist){
        name_base_dist<-names(data_base_dist)[qq]
        data_base_dist[,name_base_dist]<-mean(data_base_dist1[,name_base_dist],na.rm=T)
    }
    
    result_dist<-as.data.frame(cbind(dt<-rep(0,len_date_dist),price.pre<-rep(0,len_date_dist),r2<-rep(0,len_date_dist),FBPI<-rep(0,len_date_dist)))
    names(result_dist)<-c('DT','PRE','AdjR2','FBPI')
    
    for(jj in 1:len_date_dist){
        fit_dist<-fit.vec_dist[[pp]][[jj]]
        price.pre_dist<-predict(fit_dist,newdata=data_base_dist)
        fit_dist.sum<-summary(fit_dist)
        fit_dist.squ<-fit_dist.sum$adj.r.squared
        result_dist[jj,1]<-date_dist[jj]
        result_dist[jj,2]<-exp(price.pre_dist[[1]])
        result_dist[jj,3]<-fit_dist.squ
        result_dist[jj,4]<-result_dist[jj,2]/price.pre_base*pts_base
    }
    result.vec_dist[[pp]]<-result_dist
}

#城市与行政区最终FBPI、PRE、AdjR2
result_FBPI<-result[,c('DT','FBPI')]
names(result_FBPI)<-c('DT',city_name)

result_PRE<-result[,c('DT','PRE')]
names(result_PRE)<-c('DT',city_name)

result_AdjR2<-result[,c('DT','AdjR2')]
names(result_AdjR2)<-c('DT',city_name)

for(xx in 1:len_dist){
    dist_cd<-dist.name[xx]
    result.vec_dist.new<-result.vec_dist[names(result.vec_dist)==dist_cd][[1]]
    
    if(is.null(result.vec_dist.new)) next  #如果结果为NULL，跳出本次循环
    
    result.vec_dist.new_FBPI<-result.vec_dist.new[,c('DT','FBPI')]
    names(result.vec_dist.new_FBPI)<-c('DT',dist_cd)
    result_FBPI<-merge(result_FBPI,result.vec_dist.new_FBPI,by='DT',all=TRUE)
    
    result.vec_dist.new_PRE<-result.vec_dist.new[,c('DT','PRE')]
    names(result.vec_dist.new_PRE)<-c('DT',dist_cd)
    result_PRE<-merge(result_PRE,result.vec_dist.new_PRE,by='DT',all=TRUE)
    
    result.vec_dist.new_AdjR2<-result.vec_dist.new[,c('DT','AdjR2')]
    names(result.vec_dist.new_AdjR2)<-c('DT',dist_cd)    
    result_AdjR2<-merge(result_AdjR2,result.vec_dist.new_AdjR2,by='DT',all=TRUE)
}

#结果整理
date<-rep(0,nrow(result_FBPI))%>%as.data.frame
names(date)<-'日期'
for(yy in 1:nrow(result_FBPI)){
    int<-result_FBPI$DT[yy]%/%12
    mod<-result_FBPI$DT[yy]%%12
    if(mod==0){date$日期[yy]<-paste0(int-1,"年","12月")}else{date$日期[yy]<-paste0(int,"年",mod,"月")}
}
result_FBPI2<-cbind(date,result_FBPI)%>%select(-one_of('DT'))

if(pty==11){pty_name<-'住宅'
}else if(pty==21){pty_name<-'办公'
}else{pty_name<-'商铺'
}

# #结果保存.csv
# #保存HPI
# file_name1<-paste0(city_name,'及其行政区',pty_name,'出租')
# file.save1<-paste0("~/mnt/HPI/",file_name1,"价格指数.csv")
# write.csv(result_FBPI2,file.save1,fileEncoding = 'gbk',row.names = FALSE)
# 
# #保存PRE
# file_name2<-paste0(city_name,'及其行政区',pty_name,'出租')
# file.save2<-paste0("~/mnt/PRE/",file_name2,"pre.csv")
# write.csv(result_PRE,file.save2,fileEncoding = 'gbk',row.names = FALSE)
# 
# #保存AdjR2
# file_name3<-paste0(city_name,'及其行政区',pty_name,'出租')
# file.save3<-paste0("~/mnt/AdjR2/",file_name3,"R2.csv")
# write.csv(result_AdjR2,file.save3,fileEncoding = 'gbk',row.names = FALSE)

#结果保存.xlsx
library(xlsx)
# write.xlsx(result_FBPI2,"/home/yangsj/mnt/aaa.xlsx","Sheet1",row.names = FALSE)
file_name1<-paste0(city_name,'及其行政区',pty_name,'出租')
file.save1<-paste0("/home/yangsj/mnt/HPI/",file_name1,"价格指数.xlsx")
write.xlsx(result_FBPI2,file.save1,"Sheet1",row.names = FALSE)

#保存PRE
file_name2<-paste0(city_name,'及其行政区',pty_name,'出租')
file.save2<-paste0("/home/yangsj/mnt/PRE/",file_name2,"pre.xlsx")
write.xlsx(result_PRE,file.save2,"Sheet1",row.names = FALSE)

#保存AdjR2
file_name3<-paste0(city_name,'及其行政区',pty_name,'出租')
file.save3<-paste0("/home/yangsj/mnt/AdjR2/",file_name3,"R2.xlsx")
write.xlsx(result_AdjR2,file.save3,"Sheet1",row.names = FALSE)

